Load libraries

library(readxl)
library(tidyverse)
library(jyluMisc)
library(DrugScreenExplorer)
library(ComplexHeatmap)
library(circlize)
library(gridExtra)
library(parallel)
library(edgeR)
library(scran)
library(scater)
library(knitr)

Define variables

opt <- list()
opt$drugscreen <- "data/submission/drugScreens_pseudo.RDS"
opt$druganno <- "misc/drugList_suppl.xlsx"
opt$western <- "data/western.csv"
opt$scebulk <- "data/submission/RNA_complete.RDS"
opt$plot <- "plots/SFig6_8/"

## Set theme
lgd <-  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), 
        text = element_text(size = 15),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA), 
        legend.background = element_rect(fill='transparent',colour = NA), 
        strip.background = element_blank(), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

Load drug screen data and annotations

## Load the drug screen data
drugList <- readRDS(opt$drugscreen)
drugAnno <- read_excel(opt$druganno) %>%
  dplyr::select(-Supplier, -Screen) %>% unique()

## Load the western data
western <- read.csv(opt$western)

## RNA-Seq data
sce <- readRDS(opt$scebulk)
sce
## class: SingleCellExperiment 
## dim: 36601 117 
## metadata(0):
## assays(2): counts logcounts
## rownames(36601): ENSG00000000003 ENSG00000000005 ... ENSG00000288459
##   ENSG00000288460
## rowData names(3): ID Symbol Type
## colnames(117): H358_BafilomycinA1 H358_Birinapant ... H524_Selinexor
##   H524_Venetoclax
## colData names(5): uniqueBarcode patientID Diagnosis_simple Treatment
##   label
## reducedDimNames(2): PCA UMAP
## mainExpName: NULL
## altExpNames(0):
colnames(sce) <- sce$uniqueBarcode
## Define drug combinations
drugComb <- c("Birinapant|Necrostatin-1 (25)", "Birinapant|QVD-Oph (25)", "Birinapant|Necrostatin-1 (25)|QVD-Oph (25)", 
              "Birinapant|Necrostatin-1 (12.5)", "Birinapant|QVD-Oph (12.5)", "Birinapant|Necrostatin-1 (12.5)|QVD-Oph (12.5)",
              "GDC-0152|Necrostatin-1 (25)", "GDC-0152|QVD-Oph (25)", "GDC-0152|Necrostatin-1 (25)|QVD-Oph (25)", "GDC-0152|Necrostatin-1 (12.5)", 
              "GDC-0152|QVD-Oph (12.5)", "GDC-0152|Necrostatin-1 (12.5)|QVD-Oph (12.5)", "Birinapant|Ipatasertib (2)", 
              "Birinapant|Ruxolitinib (2)", "Birinapant|Dacinostat (0.04)", "Birinapant|Venetoclax (0.04)",
              "Birinapant|Bafilomycin_A1 (0.08)", "Birinapant|NSA (0.8)", "Birinapant|NSA (2)", 
              "Birinapant|NSA (5)", "Birinapant|NSA (12.5)", "Birinapant|NSA (0.8)|QVD-Oph (25)", 
              "Birinapant|NSA (2)|QVD-Oph (25)", "Birinapant|NSA (5)|QVD-Oph (25)", "Birinapant|NSA (12.5)|QVD-Oph (25)", 
              "Birinapant + Ipatasertib 2µM", "Birinpant + Ruxolitinib 2µM", "Birinapant + Bafilomycin A1 0.08µM", 
              "Birinapant + NSA 0.8µM",
              "Birinapant + NSA 0.8µM + QVD-Oph 25µM", "Birinapant + NSA 2µM",  "Birinapant + NSA 5µM", 
              "Birinapant + NSA 2µM + QVD-Oph 25µM", "Birinapant + NSA 5μM + QVD-Oph 25μM", 
              "DMSO", "empty", "Necrostatin-1|QVD-Oph")

SFig 6

Compute correlation of IAP-inhibitor response and protein levels

## Define proteins of interest
proteins <- colnames(western)[4:ncol(western)]

western.long <- western %>% pivot_longer(cols = c(proteins),
                         names_to = "protein", values_to = "expr") %>%
  mutate(protein = str_replace(protein, pattern = "\\.", replacement = "-")) %>%
  mutate(protein = ifelse(protein == "BCL-2", "Bcl-2", protein)) %>%
  filter(!protein %in% c("MLKL", "MLKL_upperBand")) %>%
  mutate(protein = ifelse(protein == "MLKL_lowerBand", "MLKL", protein),
         protein = ifelse(protein == "Caspase8", "Caspase-8", protein))

Scatter plots correlation of IAP-inhibitor response and protein levels

t_lymphoma <- drugList[["ScreenE"]]

## Convert proteins
proteins <- colnames(western)[4:ncol(western)]

western.long <- western %>% pivot_longer(cols = c(proteins),
                         names_to = "protein", values_to = "expr") %>%
  mutate(protein = str_replace(protein, pattern = "\\.", replacement = "-")) %>%
  mutate(protein = ifelse(protein == "BCL-2", "Bcl-2", protein)) %>%
  filter(!protein %in% c("MLKL", "MLKL_upperBand")) %>%
  mutate(protein = ifelse(protein == "MLKL_lowerBand", "MLKL", protein),
         protein = ifelse(protein == "Caspase8", "Caspase-8", protein))

## Summarise the effect over concentrations
screenDat <- t_lymphoma %>% filter(screen == "T_lymphoma") %>%
  group_by(patientID, diagnosis, screen, name) %>%
  mutate(viab.auc = mean(normVal.cor, na.rm = TRUE)) %>% ungroup()

## Get viab.auc values for birinapant and gdc-0152
drugRed <- screenDat %>% filter(diagnosis == "T-PLL", name %in% c("Birinapant", "GDC-0152")) %>%
  select(patientID, diagnosis, name, viab.auc) %>% unique()

#### Plot scatter plots for individual pairs
drugPlot <- "Birinapant"
proteinPlot <- c("RIPK3", "p53", "Bax", "Bcl-2", "cIAP-2", "cIAP-1", "Myc", "FLIPL")
  
pList <- lapply(proteinPlot, function(proteinPlot) {
  p <- drugRed %>% left_join(western.long, by = c("patientID", "diagnosis")) %>%
    filter(protein == proteinPlot, name == drugPlot) %>%
    ggplot(aes(x = viab.auc, y = expr, fill = diagnosis)) +
    ylab(paste0("Norm. ", proteinPlot, " Expression")) +
    geom_smooth(method = "lm", colour = "black", fill = "lightgrey") +
    geom_point(size = 4, shape = 21) + xlab(paste0("Avg. Viability ", drugPlot)) +
    theme_classic() +
    theme(text = element_text(size = 17.5), 
          axis.text = element_text(size = 16),
          legend.title = element_blank(),
          axis.line = element_line(linewidth = 1),
          legend.key = element_blank(),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "transparent",colour = NA),
          plot.background = element_rect(fill = "transparent",colour = NA),
          legend.background = element_rect(fill='transparent')) +
    scale_fill_manual(values = c("T-PLL" = "#64B5F6", "T-NHL" = "#A5D6A7"))
})
names(pList) <- proteinPlot
grid.arrange(grobs = pList)

lapply(proteinPlot, function(x) {
  p <- pList[[x]]
  ggsave(plot = p, filename = paste0(opt$plot, "corr_", drugPlot, "_", x, ".png"), 
         height = 4.25, width = 5.5)
})
## [[1]]
## [1] "plots/SFig6_8/corr_Birinapant_RIPK3.png"
## 
## [[2]]
## [1] "plots/SFig6_8/corr_Birinapant_p53.png"
## 
## [[3]]
## [1] "plots/SFig6_8/corr_Birinapant_Bax.png"
## 
## [[4]]
## [1] "plots/SFig6_8/corr_Birinapant_Bcl-2.png"
## 
## [[5]]
## [1] "plots/SFig6_8/corr_Birinapant_cIAP-2.png"
## 
## [[6]]
## [1] "plots/SFig6_8/corr_Birinapant_cIAP-1.png"
## 
## [[7]]
## [1] "plots/SFig6_8/corr_Birinapant_Myc.png"
## 
## [[8]]
## [1] "plots/SFig6_8/corr_Birinapant_FLIPL.png"
## Report estimates and adjusted p-values
#corTab %>% filter(name == drugPlot, protein %in% proteinPlot)

Plot RNA data

t_lymphoma <- drugList[["ScreenE"]]

## Summarise the effect over concentrations
screenDat <- t_lymphoma %>% filter(screen == "T_lymphoma") %>%
  group_by(patientID, diagnosis, screen, name) %>%
  mutate(viab.auc = mean(normVal.cor, na.rm = TRUE)) %>% ungroup()

## Select genes
geneSub <- c("BCL2", "TP53", "MLKL", "CFLAR", "BIRC2", "BIRC3", "TRAF2", "FADD", "RIPK1",
             "RIPK3", "BAX", "TNFRSF1A", "TNFRSF1B", "MYC",
             "XIAP", "CASP8", "BBC3", "BCL2A1", "BCL2L1", "MCL1", "BIK", "BID", "BAD")
colnames(sce) <- sce$uniqueBarcode

keep <- filterByExpr(assays(sce)[["counts"]], group = sce$patientID)
#keep <- which(rowData(sce)$Symbol == "ABCB1")
keep <- rowData(sce) %>% data.frame() %>%
  filter(Symbol %in% geneSub) %>% pull(ID)

dmso.mat <- logcounts(sce[keep, sce$Treatment == "DMSO"])
dmso.mat <- dmso.mat %>% data.frame() %>% rownames_to_column("gene")
colnames(dmso.mat) <- str_replace(string = colnames(dmso.mat), pattern = "_DMSO", replacement = "")

##
plotMat <- dmso.mat %>%
  left_join(select(data.frame(rowData(sce)), ID, Symbol), c("gene" = "ID")) %>%
  select(-gene) %>%
  column_to_rownames("Symbol")

plotMat <- t(scale(t(plotMat)))

## Make annotation
biri.viab <- screenDat %>% 
  filter(name == "Birinapant", patientID %in% colnames(plotMat)) %>%
  select(patientID, name, viab.auc) %>%
  unique() %>%
  pivot_wider(names_from = "patientID", values_from = "viab.auc") %>%
  column_to_rownames("name")

colAnno <- as.data.frame(matrix(NA, nrow = 1, ncol = length(colnames(plotMat))))
colnames(colAnno) <- colnames(plotMat)
name.match <- match(colnames(colAnno), colnames(biri.viab))

colAnno[, name.match[!is.na(name.match)]] <- biri.viab[, name.match[!is.na(name.match)]]

colAnno <- colAnno %>% t() %>% data.frame()
colnames(colAnno) <- "Avg.Viab."

colfun <- colorRamp2(c(min(colAnno[, 1], na.rm = TRUE), max(colAnno[, 1], na.rm = TRUE)), c("white", "darkgreen"))

pHeat <- Heatmap(plotMat, clustering_method_rows = "ward.D2", 
        clustering_method_columns = "ward.D2",
        col=colorRamp2(c(-2, 0, 2), colors = c("#1976D2", "white", "#F44336")), border = TRUE, 
        name = "Scaled Expr.", #heatmap_legend_param = list(title_gp = gpar(fontsize = 11.5, fontface="bold")),
        column_title = paste0("RNA Expression"), 
        rect_gp = gpar(col = "black", lwd = 0.1),
        column_title_gp = gpar(fontface = "plain", fontsize = 15),
        bottom_annotation = HeatmapAnnotation(df = colAnno, na_col = "grey80", 
                                           border = TRUE,
                                           annotation_label = gt_render(""), # this way the label is set to zero
                                           #gp = gpar(col = "black"), 
                                           col = list(Avg.Viab. = colfun), 
                                           which = "column")

        #col=colorRamp2(c(0, max(western.long$expr)), colors = c("grey90", "#F44336"))
        )
## Loading required namespace: gridtext
pHeat

png(file=paste0(opt$plot, "rna_heatmap.png"),
    height = 15, width = 14, res = 600, units = "cm", bg = "transparent")
draw(pHeat, background = "transparent")
dev.off()
## png 
##   2

Correlation with RNA-Seq data

#### Define genes to subset
geneSub <- c("BCL2", "TP53", "MLKL", "CFLAR", "BIRC2", "BIRC3", "TRAF2", "FADD", "RIPK1",
                       "RIPK3", "BAX", "TNFRSF1A", "TNFRSF1B", "MYC",
                       "XIAP", "CASP8", "BBC3", "BCL2A1", "BCL2L1", "MCL1", "BIK", "BID", "BAD")

screenData <- drugList[["ScreenE"]]

#keep <- filterByExpr(counts(sce), group = sce$patientID)
#keep <- which(rowData(sce)$Symbol == "ABCB1")
keep <- rowData(sce) %>% data.frame() %>%
  filter(Symbol %in% geneSub) %>% pull(ID)

dmso.mat <- logcounts(sce[keep, sce$Treatment == "DMSO"])
dmso.mat <- dmso.mat %>% data.frame() %>% rownames_to_column("gene")
colnames(dmso.mat) <- str_replace(string = colnames(dmso.mat), pattern = "_DMSO", replacement = "")

## Calculate correlation between DMSO expression and drug response
avg.df <- screenData %>% filter(!name %in% drugComb) %>%
  filter(name %in% c("Birinapant", "GDC-0152")) %>%
  group_by(patientID, name) %>%
  summarise(viab.auc = mean(normVal.cor, na.rm = TRUE)) %>%
  ungroup()
## `summarise()` has grouped output by 'patientID'. You can override using the
## `.groups` argument.
corTab <- dmso.mat %>%
    pivot_longer(cols = c(unique(sce$patientID)), 
                 names_to = "patientID", values_to = "expr") %>%
    left_join(avg.df, by = c("patientID" = "patientID")) %>%
    group_by(gene, name) %>%
    nest() %>%
    mutate(m = map(data, ~cor.test(~viab.auc + expr, .))) %>%
    mutate(res = map(m, broom::tidy)) %>%
    unnest(res) %>% ungroup() %>%
    select(name, gene, estimate, p.value) %>%
    arrange(p.value) %>%
    mutate(p.adj = p.adjust(p.value, method = "BH"))
## Warning in left_join(., avg.df, by = c(patientID = "patientID")): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 21 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
pTab <- corTab %>% left_join(select(data.frame(rowData(sce)), ID, Symbol), c("gene" = "ID")) %>%
  filter(name == "Birinapant") %>%
  filter(Symbol %in% c("BCL2", "TP53", "MLKL", "CFLAR", "BIRC2", "BIRC3", "TRAF2", "FADD", "RIPK1",
                       "RIPK3", "BAX", "TNFRSF1A", "TNFRSF1B", "MYC",
                       "XIAP", "CASP8", "BBC3", "BCL2A1", "BCL2L1", "MCL1", "BIK", "BID", "BAD", "ATM", "GAPDH")) %>%
  mutate(p.adj = p.adjust(p.value, method = "BH"))

kable(pTab)
name gene estimate p.value p.adj Symbol
Birinapant ENSG00000140379 -0.7841790 0.0072479 0.1331996 BCL2A1
Birinapant ENSG00000100290 0.7550157 0.0115826 0.1331996 BIK
Birinapant ENSG00000023445 -0.6525774 0.0408221 0.2575190 BIRC3
Birinapant ENSG00000028137 -0.6406411 0.0459673 0.2575190 TNFRSF1B
Birinapant ENSG00000003402 -0.5928923 0.0708457 0.2575190 CFLAR
Birinapant ENSG00000127191 -0.5818474 0.0776354 0.2575190 TRAF2
Birinapant ENSG00000168404 -0.5669096 0.0874680 0.2575190 MLKL
Birinapant ENSG00000137275 -0.5518125 0.0981822 0.2575190 RIPK1
Birinapant ENSG00000171791 0.5483275 0.1007683 0.2575190 BCL2
Birinapant ENSG00000067182 -0.5186120 0.1245737 0.2865195 TNFRSF1A
Birinapant ENSG00000105327 0.4796820 0.1606289 0.3358604 BBC3
Birinapant ENSG00000110330 0.4274645 0.2178667 0.4175779 BIRC2
Birinapant ENSG00000002330 0.3548211 0.3143869 0.5562230 BAD
Birinapant ENSG00000143384 -0.3348527 0.3442632 0.5655753 MCL1
Birinapant ENSG00000171552 -0.2684472 0.4532913 0.6950467 BCL2L1
Birinapant ENSG00000136997 0.2507658 0.4846626 0.6967024 MYC
Birinapant ENSG00000015475 -0.2246053 0.5327208 0.7207399 BID
Birinapant ENSG00000087088 0.1938491 0.5915335 0.7281644 BAX
Birinapant ENSG00000064012 -0.1887400 0.6015271 0.7281644 CASP8
Birinapant ENSG00000129465 0.1514256 0.6762478 0.7776850 RIPK3
Birinapant ENSG00000168040 0.1200489 0.7411451 0.8117303 FADD
Birinapant ENSG00000101966 -0.0699194 0.8477968 0.8626508 XIAP
Birinapant ENSG00000141510 -0.0630381 0.8626508 0.8626508 TP53
geneSub <- c("BCL2", "TP53", "MLKL", "CFLAR", "BIRC2", "BIRC3", "TRAF2", "FADD", "RIPK1",
                       "RIPK3", "BAX", "TNFRSF1A", "TNFRSF1B", "MYC",
                       "XIAP", "CASP8", "BBC3", "BCL2A1", "BCL2L1", "MCL1", "BIK", "BID", "BAD")

## Get into matrix format for heatmap
corMat <- corTab %>% left_join(select(data.frame(rowData(sce)), ID, Symbol), c("gene" = "ID")) %>% 
  filter(name %in% c("Birinapant", "GDC-0152"), Symbol %in% geneSub) %>% 
  mutate(p.value = -log10(p.value)) %>%
  select(-p.adj, -p.value, -gene) %>% 
  pivot_wider(names_from = "Symbol", values_from = "estimate") %>%
  column_to_rownames("name")

## Define colour scheme
colorList <- c(colorRampPalette(c("lightseagreen", "white"))(20),
               colorRampPalette(c("white"))(10),
               colorRampPalette(c("white","deeppink"))(20))

p4 <- pheatmap(t(corMat), color = colorList, 
         breaks = seq(-1,1,length.out = 50), #main = "Correlation Gene Expression and Drug Response",
                  column_title = "Corr. Gene Expression \n and Avg. Viab.",
         column_title_gp = gpar(fontface = "plain", fontsize = 15),
         border_color = "grey50",
         clustering_method = "ward.D2", name = "R^2")
p4

png(file=paste0(opt$plot, "Corr_rna_drug.png"), 
    height = 15, width = 7.5, res = 600, units = "cm", bg = "transparent")
draw(p4, background = "transparent")
dev.off()
## png 
##   2

Show scatter plots

## Show scatter plot
pList <- lapply(c("Birinapant"), function(x) {
  avg.df <- avg.df %>% filter(name == x) 
  
  pList <- lapply(c("BCL2", "BIRC2", "BIRC3", "CFLAR", "BAX", "MYC"), function(y) {
    p <- dmso.mat %>% left_join(select(data.frame(rowData(sce)), ID, Symbol), c("gene" = "ID")) %>% 
      filter(Symbol == y) %>%
      pivot_longer(cols = c(unique(sce$patientID)), 
                   names_to = "patientID", values_to = "expr") %>%
      left_join(avg.df, by = c("patientID" = "patientID")) %>%
      mutate(Diagnosis = ifelse(patientID == "H501", "T-LGL", "T-PLL")) %>%
      ggplot(aes(x = viab.auc, y = expr)) +
      geom_smooth(method = "lm", colour = "black", fill = "lightgrey") +
      geom_point(size = 4, shape = 21, aes(fill = Diagnosis)) + xlab(paste0("Avg. Viability ", x)) +
      ylab(paste0("Log Norm. Expr ", y)) +
      theme_classic() +
      theme(text = element_text(size = 17.5), 
            axis.text = element_text(size = 16),
            legend.title = element_blank(),
            axis.line = element_line(linewidth = 1),
            legend.key = element_blank(),
            panel.grid.major = element_blank(), 
            panel.grid.minor = element_blank(),
            panel.background = element_rect(fill = "transparent",colour = NA),
            plot.background = element_rect(fill = "transparent",colour = NA),
            legend.background = element_rect(fill='transparent')) +
      scale_fill_manual(values = c("T-PLL" = "#64B5F6", "T-LGL" = "#A5D6A7"))
    
    ggsave(plot = p, filename = paste0(opt$plot, "corr_rna_", x, "_", y, ".png"), 
           height = 4.25, width = 5.5)
    p
  })
  grid.arrange(grobs = pList)
})

SFig7

Overlay expected and observed viability - per patient

drug.df <- drugList[["CombiTecan"]]

drug.df <- drug.df %>% 
  separate(name, into = c("drug1", "drug2", "drug3"), sep = "[|]", remove = FALSE) %>%
  separate(drug2, into = c("drug2", "conc2"), sep = " ") %>%
  separate(drug3, into = c("drug3", "conc3"), sep = " ") %>%
  mutate(conc2 = str_replace(string = conc2, pattern = "uM", replacement = ""), 
         conc2 = as.numeric(conc2))

drug1Name <- c("Birinapant")
drug2List <- drug.df %>% filter(!drug1 %in% c("DMSO", "PBS", "empty"), !is.na(drug2), is.na(drug3)) %>% pull(drug2) %>%
  unique()

## Calculate the combination index for every drug
combiTab <- lapply(drug2List, function(x) {
  lapply(unique(drug.df$patientID), function(y) {
    drug.df <- drug.df %>% filter(patientID == y)
    
    message(x)
    # dose-response of drug1 (single agent)
    viabTab1 <- drug.df %>% filter(drug1 == drug1Name, is.na(drug2)) %>% dplyr::rename(viab1 = normVal)
  
    # dose-response of drug1 in combination with drug2
    viabTab2 <-  drug.df %>% filter(drug1 == drug1Name, drug2 %in% x, is.na(drug3)) %>% dplyr::rename(viab2 = normVal) 
  
    # get the viability table of the drug with only 1 concentration (base drug)
    viabSingle <- drug.df %>% filter(drug1 == x, concentration == unique(viabTab2$conc2), is.na(drug2)) %>%
      distinct(patientID, normVal) %>% dplyr::rename(viab3 = normVal)

    testTab <- viabTab2 %>%
          left_join(select(viabTab1, patientID, viab1, concentration)) %>%
          left_join(viabSingle) %>%
          mutate(expViab = viab1*viab3) %>%
          mutate(CI = viab2/expViab, meanViab2 = mean(viab2)) %>%
          filter(!is.na(viab1), !is.na(viab2),!is.na(viab3)) %>%
    dplyr::rename(obsViab = viab2)
  }) %>% bind_rows()
}) %>% bind_rows()

lapply(c("Necrostatin-1", "NSA", "GSK-872", "QVD-Oph"), function(x) {
  concInter <- drug.df %>% filter(name == "Birinapant") %>% pull(concentration) %>% unique()
  
  pList <- lapply(unique(combiTab$patientID), function(y) {
    combiTab <- combiTab %>% filter(patientID == y) %>% 
      filter(drug2 %in% x)
    
    if(x == "QVD-Oph")  {lim.plot <- 1.6} else {lim.plot <- 1.25} 
  
    p <- combiTab %>%
      dplyr::rename(Expected = expViab, Observed = obsViab, Single = viab1) %>%
      pivot_longer(cols = c(Expected, Observed, Single), names_to = "viab",
                   values_to = "values") %>%
      mutate(viab = factor(viab, levels = c("Observed", "Expected", "Single")),
             drug2 = factor(drug2, levels = x), 
             patientID = paste0(patientID, " (", diagnosis, ")")) %>%
      arrange(desc(viab)) %>%
      #filter(viab %in% c("Observed", "Expected")) %>%
      ggplot(aes(x = concentration, y = values, group = viab)) +
      geom_line(linewidth = 0.5, aes(colour = viab)) +
      geom_point(size = 3.5, shape = 21, aes(fill = viab)) +
      geom_hline(yintercept = 1, linetype = "dashed", colour = "black", linewidth = 0.75) +
      geom_hline(yintercept = unique(combiTab$viab3), linetype = "dashed", colour = "darkorange", linewidth = 0.75) +
      scale_x_log10(labels = concInter, breaks = concInter) +
      scale_y_continuous(limits = c(0, lim.plot), breaks = c(0, 0.5, 1.0, 1.5)) +
      xlab("Birinapant (µM)") + ylab("Viability") + #ggtitle(x) +
      theme_classic() +
      theme(text = element_text(size = 22.5),
            plot.title = element_text(hjust = 0.5, size = 22.5),
            panel.background = element_rect(fill = "transparent",colour = NA),
            plot.background = element_rect(fill = "transparent",colour = NA),
            legend.background = element_rect(fill='transparent'),
            strip.background = element_blank(),
            #strip.text.x = element_text(size = 12.5),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.title = element_blank(),
            axis.line = element_line(linewidth = 1),
            #axis.text.y = element_text(size = 12.5),
            #axis.text.x = element_text(angle = 45, hjust=1, vjust = 1),
            legend.position = "none"
      ) +
      #scale_color_manual(values = c("Single" = "black", "Observed" = "#FFA000", "Expected" = "#4DB6AC"))
      scale_color_manual(values = c("Observed" = "#F44336", "Expected" = "#1976D2", "Single" = "black")) +
      scale_fill_manual(values = c("Observed" = "#F44336", "Expected" = "#1976D2", "Single" = "black")) +
      facet_wrap(. ~ patientID)

    # ggsave(plot = p, filename = paste0(opt$plot, x, "_combi_curve.png"),
    #        height = 4, width = 4.75)

    p
  })
  ggsave(plot = grid.arrange(grobs = pList, name = x), filename = paste0(opt$plot, x, "_combi_curve.png"),
         height = 35, width = 30, units = "cm")
})

## [[1]]
## [1] "plots/SFig6_8/Necrostatin-1_combi_curve.png"
## 
## [[2]]
## [1] "plots/SFig6_8/NSA_combi_curve.png"
## 
## [[3]]
## [1] "plots/SFig6_8/GSK-872_combi_curve.png"
## 
## [[4]]
## [1] "plots/SFig6_8/QVD-Oph_combi_curve.png"

Scatter plot birinapant + QVD-Oph vs birinapant + necrostatin-1

screenData <- drugList[["ScreenE"]]

screenData <- screenData %>% select(screen, wellID, value, patientID, sampleID, diagnosis, name,
                                    concentration, wellType, normVal, normVal.cor, concIndex) %>%
  mutate(name = str_replace(name, "Bafilomycin A1", "Bafilomycin_A1")) %>%
  filter(screen == "T_lymphoma") %>%
  separate(name, into = c("drug1","drugB","drugC"), sep = "[|]", remove = FALSE) %>%
  separate(drugB, into=c("drug2", "conc2"), sep = " ") %>%
  separate(drugC , into = c("drug3", "conc3"), sep = " ") %>%
  mutate(conc2 = as.numeric(str_remove_all(conc2,"[() ]")),
         conc3 = as.numeric(str_remove_all(conc3, "[() ]"))) %>%
  mutate(conc2 = ifelse(drug2 == "QVD-Oph" & is.na(conc2), concentration, conc2))

comObs <- screenData %>% filter(drug1 %in% c("Birinapant", "GDC-0152"), !is.na(drug2), is.na(drug3)) %>% #only 2 drug combinations
  select(patientID, sampleID, diagnosis, name, drug1, drug2, concentration, conc2, normVal.cor) %>%
  dplyr::rename(comViab = normVal.cor)

drugObs <- screenData %>% filter(drug1 %in% c("Birinapant", "GDC-0152"), is.na(drug2)) %>%
  select(sampleID, drug1, concentration, normVal.cor) %>% 
  dplyr::rename(drugViab = normVal.cor)

baseObs <- screenData %>% filter(drug1 %in% c("Ipatasertib", "NSA", 
                                              "Bafilomycin_A1", "Ruxolitinib", 
                                              "Venetoclax", "Dacinostat", "Necrostatin-1", "QVD-Oph"), is.na(drug2)) %>%
  select(sampleID, drug1, concentration, normVal.cor) %>%
  dplyr::rename(drug2 = drug1, conc2 = concentration, baseViab = normVal.cor)

dmsoObs <- screenData %>% filter(drug1 == "DMSO") %>%
  group_by(sampleID) %>% summarise(dmsoViab = mean(normVal.cor))

## CombiConc
combiConc <- 12.5

combiTab <- left_join(comObs, drugObs, by = c("sampleID", "drug1", "concentration")) %>%
  left_join(baseObs, by = c("sampleID","drug2", "conc2")) %>%
  left_join(dmsoObs, by = "sampleID") %>%
#  mutate(expViab = drugViab*baseViab/dmsoViab) %>%
  mutate(expViab = drugViab*baseViab) %>%
  mutate(CI =  comViab/expViab) %>%
  arrange(name) %>%
  filter(conc2 %in% combiConc)

## Drug combinations necrostatin-1 and QVD-Oph
sumCombiTab <- combiTab %>% 
  filter(drug2 %in% c("Necrostatin-1", "QVD-Oph")) %>%
  mutate(diagnosis = ifelse(diagnosis %in% c("PTCL", "T-LGL", "Sezary", "AITL"), "T-NHL", diagnosis)) %>%
  mutate(diagnosis = ifelse(diagnosis %in% c("CLL", "MCL", "MZL", "DLBCL", "HCL"), "B-cell", diagnosis)) %>%
  group_by(patientID, diagnosis, drug1, drug2) %>%
  summarise(meanCI = log(mean(CI))) %>% ungroup()
                       
plotTab <- sumCombiTab %>% 
  pivot_wider(names_from = drug2, values_from = meanCI)
                       
p5 <- plotTab %>%
  filter(drug1 == "Birinapant") %>%
  dplyr::rename(Diagnosis = diagnosis) %>%
  ggplot(aes(x=`Necrostatin-1`, y=`QVD-Oph`, fill = Diagnosis, label = patientID)) +
    geom_point(size = 4, shape = 21, stroke = 0.2) +
    geom_hline(yintercept = 0) + geom_vline(xintercept = 0) +
  xlab("Log2(Avg. CI) - Necrostatin-1") +
  ylab("Log2(Avg. CI) - Q-VD-Oph") +
    ggtitle("Birinapant") +
    theme_classic() +
    theme(text = element_text(size = 17.5),
          plot.title = element_text(hjust = 0.5),
     # axis.text=element_text(size=12),
      #axis.title=element_text(size=14),
      #legend.position = "none", 
      strip.background = element_blank(),
      #strip.text = element_text(size = 12), 
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.line = element_line(linewidth = 1),
      panel.background = element_rect(fill = "transparent",colour = NA),
      plot.background = element_rect(fill = "transparent",colour = NA), 
      legend.background = element_rect(fill='transparent')) +
    scale_fill_manual(values = c("B-cell" = "#F44336", "T-PLL" = "#64B5F6", "T-NHL" = "#A5D6A7")) +
    #scale_fill_manual(values = c("B-cell" = "#F44336", "T-PLL" = "#CE93D8", "Other T-cell" = "#FDD835")) +
  guides(fill = guide_legend(override.aes = list(size=3.5)))
p5

ggsave(plot = p5, filename = paste0(opt$plot, "QVD_vs_Necro_biri.png"), 
       height = 5, width = 6.5)

Scatter plot birinapant + QVD-Oph vs birinapant + NSA

screenData <- drugList[["ScreenE"]]

screenData <- screenData %>% select(screen, wellID, value, patientID, sampleID, diagnosis, name,
                                    concentration, wellType, normVal, concIndex) %>%
  unique() %>%
  mutate(diagnosis = ifelse(diagnosis %in% c("Sezary", "T-LGL", "AITL", "PTCL"), "T-NHL", diagnosis)) %>%
  mutate(name = str_replace(name, "Bafilomycin A1", "Bafilomycin_A1")) %>%
  filter(screen == "T_lymphoma_combi") %>%
  separate(name, into = c("drug1","drugB","drugC"), sep = "[|]", remove = FALSE) %>%
  separate(drugB, into=c("drug2", "conc2"), sep = " ") %>%
  separate(drugC , into = c("drug3", "conc3"), sep = " ") %>%
  mutate(conc2 = as.numeric(str_remove_all(conc2,"[() ]")),
         conc3 = as.numeric(str_remove_all(conc3, "[() ]"))) %>%
  mutate(conc2 = ifelse(drug2 == "QVD-Oph" & is.na(conc2), concentration, conc2))
## Warning: Expected 3 pieces. Missing pieces filled with `NA` in 3192 rows [1, 2, 3, 4, 5,
## 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].
comObs <- screenData %>% filter(drug1 %in% c("Birinapant", "GDC-0152"), !is.na(drug2), is.na(drug3)) %>% #only 2 drug combinations
  select(patientID, sampleID, diagnosis, name, drug1, drug2, concentration, conc2, normVal) %>%
  dplyr::rename(comViab = normVal) %>% unique()

drugObs <- screenData %>% filter(drug1 %in% c("Birinapant", "GDC-0152"), is.na(drug2)) %>%
  select(sampleID, drug1, concentration, normVal) %>% 
  dplyr::rename(drugViab = normVal) %>% unique()

baseObs <- screenData %>% filter(drug1 %in% c("NSA", "Necrostatin-1", "QVD-Oph"), is.na(drug2)) %>%
  select(sampleID, drug1, concentration, normVal) %>%
  dplyr::rename(drug2 = drug1, conc2 = concentration, baseViab = normVal) %>% unique()

dmsoObs <- screenData %>% filter(drug1 == "DMSO") %>%
  group_by(sampleID) %>% summarise(dmsoViab = mean(normVal))

combiTabNSA <- left_join(comObs, drugObs, by = c("sampleID", "drug1", "concentration")) %>%
  left_join(baseObs, by = c("sampleID","drug2", "conc2")) %>%
  left_join(dmsoObs, by = "sampleID") %>%
#  mutate(expViab = drugViab*baseViab/dmsoViab) %>%
  mutate(expViab = drugViab*baseViab) %>%
  mutate(CI =  comViab/expViab) %>%
  arrange(name)

combiTabNSA
## # A tibble: 280 × 14
##    patientID sampleID    diagnosis name  drug1 drug2 concentration conc2 comViab
##    <chr>     <chr>       <chr>     <chr> <chr> <chr>         <dbl> <dbl>   <dbl>
##  1 H023      11PB0037_R1 CLL       Biri… Biri… NSA          10         2   0.459
##  2 H023      11PB0037_R1 CLL       Biri… Biri… NSA           2         2   0.572
##  3 H023      11PB0037_R1 CLL       Biri… Biri… NSA           0.4       2   0.678
##  4 H023      11PB0037_R1 CLL       Biri… Biri… NSA           0.08      2   0.725
##  5 H023      11PB0037_R1 CLL       Biri… Biri… NSA           0.016     2   0.865
##  6 H013      11PB0358_R1 CLL       Biri… Biri… NSA          10         2   0.718
##  7 H013      11PB0358_R1 CLL       Biri… Biri… NSA           2         2   0.867
##  8 H013      11PB0358_R1 CLL       Biri… Biri… NSA           0.4       2   0.937
##  9 H013      11PB0358_R1 CLL       Biri… Biri… NSA           0.08      2   0.964
## 10 H013      11PB0358_R1 CLL       Biri… Biri… NSA           0.016     2   0.969
## # ℹ 270 more rows
## # ℹ 5 more variables: drugViab <dbl>, baseViab <dbl>, dmsoViab <dbl>,
## #   expViab <dbl>, CI <dbl>
drug1Name <- c("Birinapant")
drug2List <- c("NSA")
#diagList <- c("T-PLL", "CLL", "MCL", "Sezary", "T-LGL", "PTCL")
diagList <- c("CLL", "MCL", "T-PLL", "T-NHL")
combiConc <- c(2)

combiTestResNSA <- lapply(drug1Name, function(y) {
  lapply(drug2List, function(x) {
    res <- lapply(diagList, function(diag) {
      message(diag)
      screenSub <- filter(screenData, diagnosis == diag, is.na(drug3))
      testTab <- lapply(combiConc, function(conc){
        viabTab1 <- filter(screenSub, drug1 == y, is.na(drug2)) %>%
          dplyr::rename(viabBackbone = normVal) # get viability of backbone drug
        viabTab2 <- filter(screenSub, drug1 == y, drug2 %in% x, conc2 == conc) %>%
          dplyr::rename(viabCombi = normVal) # get viability of combination (at different concentrations)
        viabTab3 <- filter(screenSub, drug1 == x, concentration == conc, is.na(drug2)) %>%
          dplyr::rename(viabSingle = normVal) # get viability of comb, drug as single compound.
        mergeTab <- viabTab2 %>%
          left_join(select(viabTab1, sampleID, viabBackbone, concentration)) %>%
          left_join(select(viabTab3, sampleID, viabSingle))  %>%
          mutate(expViab = viabBackbone*viabSingle) %>%
          mutate(CI = viabCombi/expViab) %>%
          filter(!is.na(viabBackbone), !is.na(viabCombi),!is.na(viabSingle))
      }) %>% bind_rows()
      lapply(unique(testTab$concentration), function(concBackbone) {
        eachConcTab <- filter(testTab, concentration == concBackbone)
        sumTab <- eachConcTab %>%
          group_by(patientID, concentration) %>%
          mutate(meanExpViab = mean(expViab), meanViabCombi = mean(viabCombi)) %>% # Summarise over combi concentrations
          ungroup() %>% select(patientID, diagnosis, drug1, drug2, concentration, meanExpViab, meanViabCombi) %>%
          unique()
        prepTab <- sumTab %>%
          pivot_longer(cols = c("meanExpViab", "meanViabCombi"), values_to = "response", names_to = "type")
        p.val <- t.test(response ~ type, data = prepTab, equal.var = TRUE, paired = TRUE)$p.value
        meanTab <- eachConcTab %>%
          group_by(concentration) %>%
          mutate(meanCI = mean(CI)) %>% select(meanCI, concentration) %>% unique()
        tibble(drug1 = y, drug2 = x, diagnosis = diag, concentration = concBackbone, p.value = p.val,
               meanCI = meanTab$meanCI)
      }) %>% bind_rows()
    }) %>% bind_rows()
  }) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p.value, method = "BH"))
## CLL
## Joining with `by = join_by(sampleID, concentration)`
## Joining with `by = join_by(sampleID)`
## MCL
## Joining with `by = join_by(sampleID, concentration)`
## Joining with `by = join_by(sampleID)`
## T-PLL
## Joining with `by = join_by(sampleID, concentration)`
## Joining with `by = join_by(sampleID)`
## T-NHL
## Joining with `by = join_by(sampleID, concentration)`
## Joining with `by = join_by(sampleID)`
sumCombiTab <- rbind(combiTab, combiTabNSA) %>% filter(drug2 %in% c("NSA", "QVD-Oph")) %>%
  mutate(diagnosis = ifelse(diagnosis %in% c("PTCL", "T-LGL", "Sezary", "AITL"), "T-NHL", diagnosis)) %>%
  mutate(diagnosis = ifelse(diagnosis %in% c("CLL", "MCL", "MZL", "DLBCL"), "B-cell", diagnosis)) %>%
  group_by(patientID, diagnosis, drug1, drug2) %>%
  summarise(meanCI = log(mean(CI))) %>% ungroup()
## `summarise()` has grouped output by 'patientID', 'diagnosis', 'drug1'. You can
## override using the `.groups` argument.
plotTab <- sumCombiTab %>%
  pivot_wider(names_from = drug2, values_from = meanCI)

p4 <- plotTab %>%
  filter(drug1 == "Birinapant") %>%
  dplyr::rename(Diagnosis = diagnosis) %>%
  ggplot(aes(x=`NSA`, y=`QVD-Oph`, fill = Diagnosis, label = patientID)) +
  geom_point(size = 4, shape = 21, stroke = 0.2) +
  geom_hline(yintercept = 0) + geom_vline(xintercept = 0) +
  xlab("Log2(Avg. CI) - NSA") +
  ylab("Log2(Avg. CI) - Q-VD-Oph") +
  ggtitle("Birinapant") +
  theme_classic() +
  theme(text = element_text(size = 17.5),
        plot.title = element_text(hjust = 0.5),
        # axis.text=element_text(size=12),
        #axis.title=element_text(size=14),
        #legend.position = "none",
        strip.background = element_blank(),
        #strip.text = element_text(size = 12),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(linewidth = 1),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA),
        legend.background = element_rect(fill='transparent')) +
  scale_fill_manual(values = c("B-cell" = "#F44336", "T-PLL" = "#64B5F6", "T-NHL" = "#A5D6A7"))
p4
## Warning: Removed 33 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave(plot = p4, filename = paste0(opt$plot, "QVD_vs_NSA.png"),
       height = 5, width = 6.5)
## Warning: Removed 33 rows containing missing values or values outside the scale range
## (`geom_point()`).

SFig8

Heatmap GDC-0152 + QVD-Oph/necrostatin-1

t_lymphoma <- drugList[["ScreenE"]]

## Drug combinations (with normVal.cor)
screenData <- t_lymphoma %>%
  mutate(diagnosis = ifelse(diagnosis %in% c("Sezary", "T-LGL", "AITL", "PTCL"), "T-NHL", diagnosis))

screenData <- screenData %>% select(screen, wellID, value, patientID, sampleID, diagnosis, name,
                                    concentration, wellType, normVal, normVal.cor, concIndex) %>%
  mutate(name = str_replace(name, "Bafilomycin A1", "Bafilomycin_A1")) %>%
  filter(screen == "T_lymphoma") %>%
  separate(name, into = c("drug1","drugB","drugC"), sep = "[|]", remove = FALSE) %>%
  separate(drugB, into=c("drug2", "conc2"), sep = " ") %>%
  separate(drugC , into = c("drug3", "conc3"), sep = " ") %>%
  mutate(conc2 = as.numeric(str_remove_all(conc2,"[() ]")),
         conc3 = as.numeric(str_remove_all(conc3, "[() ]"))) %>%
  mutate(conc2 = ifelse(drug2 == "QVD-Oph" & is.na(conc2), concentration, conc2))


drug1Name <- c("GDC-0152")
drug2List <- c("Necrostatin-1", "QVD-Oph")
#diagList <- c("T-PLL", "MCL", "CLL", "T-LGL", "Sezary", "PTCL") 
diagList <- c("CLL", "MCL", "T-PLL", "T-NHL") 
combiConc <- c(12.5#, 25
               )

## paired t-test (mean of drug 2 conc)
combiTestRes <- lapply(drug1Name, function(y) {  
  lapply(drug2List, function(x) {
    res <- lapply(diagList, function(diag) {
      screenSub <- filter(screenData, diagnosis == diag, is.na(drug3))
      testTab <- lapply(combiConc, function(conc){
        viabTab1 <- filter(screenSub, drug1 == y, is.na(drug2)) %>%
          dplyr::rename(viabBackbone = normVal.cor) # get viability of backbone drug
        viabTab2 <- filter(screenSub, drug1 == y, drug2 %in% x, conc2 == conc) %>%
          dplyr::rename(viabCombi = normVal.cor) # get viability of combination (at different concentrations)
        viabTab3 <- filter(screenSub, drug1 == x, concentration == conc, is.na(drug2)) %>%
          dplyr::rename(viabSingle = normVal.cor) # get viability of comb, drug as single compound.
        mergeTab <- viabTab2 %>%
          left_join(select(viabTab1, sampleID, viabBackbone, concentration)) %>%
          left_join(select(viabTab3, sampleID, viabSingle))  %>%
          mutate(expViab = viabBackbone*viabSingle) %>%
          mutate(CI = viabCombi/expViab) %>%
          filter(!is.na(viabBackbone), !is.na(viabCombi),!is.na(viabSingle))
      }) %>% bind_rows()
      lapply(unique(testTab$concentration), function(concBackbone) {
        eachConcTab <- filter(testTab, concentration == concBackbone)
        sumTab <- eachConcTab %>%
          group_by(patientID, concentration) %>%
          mutate(meanExpViab = mean(expViab), meanViabCombi = mean(viabCombi)) %>% # Summarise over combi concentrations
          ungroup() %>% select(patientID, diagnosis, drug1, drug2, concentration, meanExpViab, meanViabCombi) %>% 
          unique()
        prepTab <- sumTab %>% 
          pivot_longer(cols = c("meanExpViab", "meanViabCombi"), values_to = "response", names_to = "type")
        p.val <- t.test(response ~ type, data = prepTab, equal.var = TRUE, paired = TRUE)$p.value
        meanTab <- eachConcTab %>%
          group_by(concentration) %>%
          mutate(meanCI = mean(CI)) %>% select(meanCI, concentration) %>% unique()
        tibble(drug1 = y, drug2 = x, diagnosis = diag, concentration = concBackbone, p.value = p.val, 
               meanCI = meanTab$meanCI)
      }) %>% bind_rows()
    }) %>% bind_rows()
  }) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p.value, method = "BH"))

plotTab <- combiTestRes %>%
  mutate(pSign = -log10(p.adj)*sign(log(meanCI)), 
         sigSign = ifelse(p.adj < 0.05, "*","")) 

p4 <- plotTab %>%
  mutate(concentration = round(concentration, 5), 
         drug2 = ifelse(drug2 == "QVD-Oph", "Q-VD-Oph", drug2),
         #drug2 = factor(drug2, levels = c("Q-VD-Oph", "Necrostatin-1")), 
         diagnosis = factor(diagnosis, levels = diagList),
         drug2 = paste0(drug1, " + ", drug2)) %>%
  ggplot(aes(x = diagnosis, y = factor(concentration), fill = pSign)) +
  geom_tile(colour = "black", size = 0.05) +
  xlab("") +
  ylab("Concentration (µM)") +
  facet_wrap(. ~ drug2, scale = "free") +
  scale_fill_gradient2(high = "#F44336", mid = "white", low = "#1976D2", midpoint = 0,
                       name = "-Log10(p.adj) \n with Direction") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, hjust=1, vjust = 0.5, size = 11.5),
        axis.text.y = element_text(size = 11.5), 
        axis.title = element_text(size=13.5),
        strip.background = element_blank(), 
        strip.text.x = element_text(size = 11.5), 
        line = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA),
        legend.title = element_text(size = 12.5),
        legend.text = element_text(size = 11.5),
        legend.background = element_rect(fill='transparent'), 
        #panel.border = element_rect(colour = "black", fill=NA, size=1)
        panel.border = element_rect(colour = "black", fill=NA, size=1))
p4

ggsave(plot = p4, paste0(opt$plot, "CI_", drug1Name, "_heatmap.png"), height = 4, width = 8.25)

Calculating the combination effect

t_lymphoma <- drugList[["ScreenE"]]

## Drug combinations (with normVal.cor)
screenData <- t_lymphoma %>%
  mutate(diagnosis = ifelse(diagnosis %in% c("Sezary", "T-LGL", "AITL", "PTCL"), "T-NHL", diagnosis))
screenData <- screenData %>% select(screen, wellID, value, patientID, sampleID, diagnosis, name,
                                    concentration, wellType, normVal, normVal.cor, concIndex) %>%
  mutate(name = str_replace(name, "Bafilomycin A1", "Bafilomycin_A1")) %>%
  filter(screen == "T_lymphoma") %>%
  separate(name, into = c("drug1","drugB","drugC"), sep = "[|]", remove = FALSE) %>%
  separate(drugB, into=c("drug2", "conc2"), sep = " ") %>%
  separate(drugC , into = c("drug3", "conc3"), sep = " ") %>%
  mutate(conc2 = as.numeric(str_remove_all(conc2,"[() ]")),
         conc3 = as.numeric(str_remove_all(conc3, "[() ]"))) %>%
  mutate(conc2 = ifelse(drug2 == "QVD-Oph" & is.na(conc2), concentration, conc2))

comObs <- screenData %>% filter(drug1 %in% c("Birinapant", "GDC-0152"), !is.na(drug2), is.na(drug3)) %>% #only 2 drug combinations
  select(patientID, sampleID, diagnosis, name, drug1, drug2, concentration, conc2, normVal.cor) %>%
  dplyr::rename(comViab = normVal.cor)

drugObs <- screenData %>% filter(drug1 %in% c("Birinapant", "GDC-0152"), is.na(drug2)) %>%
  select(sampleID, drug1, concentration, normVal.cor) %>% 
  dplyr::rename(drugViab = normVal.cor)

baseObs <- screenData %>% filter(drug1 %in% c("Ipatasertib", "NSA", 
                                              "Bafilomycin_A1", "Ruxolitinib", 
                                              "Venetoclax", "Dacinostat", "Necrostatin-1", "QVD-Oph"), is.na(drug2)) %>%
  select(sampleID, drug1, concentration, normVal.cor) %>%
  dplyr::rename(drug2 = drug1, conc2 = concentration, baseViab = normVal.cor)

dmsoObs <- screenData %>% filter(drug1 == "DMSO") %>%
  group_by(sampleID) %>% summarise(dmsoViab = mean(normVal.cor))

combiConc <- 12.5

combiTab <- left_join(comObs, drugObs, by = c("sampleID", "drug1", "concentration")) %>%
  left_join(baseObs, by = c("sampleID","drug2", "conc2")) %>%
  left_join(dmsoObs, by = "sampleID") %>%
#  mutate(expViab = drugViab*baseViab/dmsoViab) %>%
  mutate(expViab = drugViab*baseViab) %>%
  mutate(CI =  comViab/expViab) %>%
  arrange(name) %>%
  filter(conc2 %in% combiConc)

Scatter plot GDC-0152 + QVD-Oph vs birinapant + necrostatin-1

## Drug combinations necrostatin-1 and QVD-Oph
sumCombiTab <- combiTab %>% filter(conc2 == combiConc) %>%
  filter(drug2 %in% c("Necrostatin-1", "QVD-Oph")) %>%
  mutate(diagnosis = ifelse(diagnosis %in% c("PTCL", "T-LGL", "Sezary", "AITL"), "T-NHL", diagnosis)) %>%
  mutate(diagnosis = ifelse(diagnosis %in% c("CLL", "MCL", "MZL", "DLBCL"), "B-cell", diagnosis)) %>%
  group_by(patientID, diagnosis, drug1, drug2) %>%
  summarise(meanCI = log(mean(CI))) %>% ungroup()
                       
plotTab <- sumCombiTab %>% 
  pivot_wider(names_from = drug2, values_from = meanCI)

p4 <- plotTab %>%
  filter(drug1 == "GDC-0152") %>%
  dplyr::rename(Diagnosis = diagnosis) %>%
  ggplot(aes(x=`Necrostatin-1`, y=`QVD-Oph`, fill = Diagnosis, label = patientID)) +
  geom_point(size = 4, shape = 21, stroke = 0.2) +
  geom_hline(yintercept = 0) + geom_vline(xintercept = 0) +
  ggtitle("GDC-0152") +
  xlab("Log2(Avg. CI) - Necrostatin-1") +
  ylab("Log2(Avg. CI) - Q-VD-Oph") +
  theme_classic() +
  theme(text = element_text(size = 17.5),
        plot.title = element_text(hjust = 0.5),
        # axis.text=element_text(size=12),
        #axis.title=element_text(size=14),
        #legend.position = "none", 
        strip.background = element_blank(),
        #strip.text = element_text(size = 12), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.line = element_line(linewidth = 0.75),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA), 
        legend.background = element_rect(fill='transparent')) +
  scale_fill_manual(values = c("B-cell" = "#F44336", "T-PLL" = "#64B5F6", "T-NHL" = "#A5D6A7")) +
  #scale_fill_manual(values = c("B-cell" = "#F44336", "T-PLL" = "#CE93D8", "Other T-cell" = "#FDD835")) +
  guides(fill = guide_legend(override.aes = list(size=3.5)))

p4

ggsave(plot = p4, paste0(opt$plot, "QVD_vs_Necro_gdc.png"), height = 5, width = 6.5)

Output per patient dose response curves - T-PLL

drug.df <- drugList[["CombiTecan"]]

drug.df <- drug.df %>% 
  separate(name, into = c("drug1", "drug2", "drug3"), sep = "[|]", remove = FALSE) %>%
  separate(drug2, into = c("drug2", "conc2"), sep = " ") %>%
  separate(drug3, into = c("drug3", "conc3"), sep = " ") %>%
  mutate(conc2 = str_replace(string = conc2, pattern = "uM", replacement = ""), 
         conc2 = as.numeric(conc2))

drug1Name <- c("Birinapant")
drug2List <- drug.df %>% filter(!drug1 %in% c("DMSO", "PBS", "empty"), !is.na(drug2), is.na(drug3)) %>% pull(drug2) %>%
  unique()

## Calculate the combination index for every drug
combiTab <- lapply(drug2List, function(x) {
  lapply(unique(drug.df$patientID), function(y) {
    drug.df <- drug.df %>% filter(patientID == y)
    
    message(x)
    # dose-response of drug1 (single agent)
    viabTab1 <- drug.df %>% filter(drug1 == drug1Name, is.na(drug2)) %>% dplyr::rename(viab1 = normVal)
  
    # dose-response of drug1 in combination with drug2
    viabTab2 <-  drug.df %>% filter(drug1 == drug1Name, drug2 %in% x, is.na(drug3)) %>% dplyr::rename(viab2 = normVal) 
  
    # get the viability table of the drug with only 1 concentration (base drug)
    viabSingle <- drug.df %>% filter(drug1 == x, concentration == unique(viabTab2$conc2), is.na(drug2)) %>%
      distinct(patientID, normVal) %>% dplyr::rename(viab3 = normVal)

    testTab <- viabTab2 %>%
          left_join(select(viabTab1, patientID, viab1, concentration)) %>%
          left_join(viabSingle) %>%
          mutate(expViab = viab1*viab3) %>%
          mutate(CI = viab2/expViab, meanViab2 = mean(viab2)) %>%
          filter(!is.na(viab1), !is.na(viab2),!is.na(viab3)) %>%
    dplyr::rename(obsViab = viab2)
  }) %>% bind_rows()
}) %>% bind_rows()

lapply(c("Ibrutinib", "Motolimod"), function(x) {
  concInter <- drug.df %>% filter(name == "Birinapant") %>% pull(concentration) %>% unique()
  
  pList <- lapply(unique(combiTab$patientID), function(y) {
    combiTab <- combiTab %>% filter(patientID == y) %>% 
      filter(drug2 %in% x)
    
    if(x == "Motolimod")  {
      lim.plot <- 1.5
    } else if(x == "QVD-Oph") {
      lim.plot <- 1.5
      } else {lim.plot <- 1.25}

    p <- combiTab %>%
      dplyr::rename(Expected = expViab, Observed = obsViab, Single = viab1) %>%
      pivot_longer(cols = c(Expected, Observed, Single), names_to = "viab",
                   values_to = "values") %>%
      mutate(viab = factor(viab, levels = c("Observed", "Expected", "Single")),
             drug2 = factor(drug2, levels = x), 
             patientID = paste0(patientID, " (", diagnosis, ")")) %>%
      arrange(desc(viab)) %>%
      #filter(viab %in% c("Observed", "Expected")) %>%
      ggplot(aes(x = concentration, y = values, group = viab)) +
      geom_line(linewidth = 0.5, aes(colour = viab)) +
      geom_point(size = 3.5, shape = 21, aes(fill = viab)) +
      geom_hline(yintercept = 1, linetype = "dashed", colour = "black", linewidth = 0.75) +
      geom_hline(yintercept = unique(combiTab$viab3), linetype = "dashed", colour = "darkorange", linewidth = 0.75) +
      scale_x_log10(labels = concInter, breaks = concInter) +
      scale_y_continuous(limits = c(0, lim.plot), breaks = c(0, 0.5, 1.0, 1.5)) +
      xlab("Birinapant (µM)") + ylab("Viability") + #ggtitle(x) +
      facet_wrap(. ~ patientID) +
      theme_classic() +
      theme(text = element_text(size = 22.5),
            plot.title = element_text(hjust = 0.5, size = 22.5),
            panel.background = element_rect(fill = "transparent",colour = NA),
            plot.background = element_rect(fill = "transparent",colour = NA),
            legend.background = element_rect(fill='transparent'),
            strip.background = element_blank(),
            #strip.text.x = element_text(size = 12.5),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.title = element_blank(),
            axis.line = element_line(linewidth = 1),
            #axis.text.y = element_text(size = 12.5),
            #axis.text.x = element_text(angle = 45, hjust=1, vjust = 1),
            #legend.position = c(0.25, 0.25
            legend.position = "none"
      ) +
      #scale_color_manual(values = c("Single" = "black", "Observed" = "#FFA000", "Expected" = "#4DB6AC"))
      scale_color_manual(values = c("Observed" = "#F44336", "Expected" = "#1976D2", "Single" = "black")) +
      scale_fill_manual(values = c("Observed" = "#F44336", "Expected" = "#1976D2", "Single" = "black"))
      #p

    #p
  })
  ggsave(plot = grid.arrange(grobs = pList, name = x), filename = paste0(opt$plot, x, "_combi_curve.png"),
         height = 35, width = 30, units = "cm")
  
  grid.arrange(grobs = pList, name = x)
})

## [[1]]
## TableGrob (4 x 3) "Ibrutinib": 12 grobs
##     z     cells      name           grob
## 1   1 (1-1,1-1) Ibrutinib gtable[layout]
## 2   2 (1-1,2-2) Ibrutinib gtable[layout]
## 3   3 (1-1,3-3) Ibrutinib gtable[layout]
## 4   4 (2-2,1-1) Ibrutinib gtable[layout]
## 5   5 (2-2,2-2) Ibrutinib gtable[layout]
## 6   6 (2-2,3-3) Ibrutinib gtable[layout]
## 7   7 (3-3,1-1) Ibrutinib gtable[layout]
## 8   8 (3-3,2-2) Ibrutinib gtable[layout]
## 9   9 (3-3,3-3) Ibrutinib gtable[layout]
## 10 10 (4-4,1-1) Ibrutinib gtable[layout]
## 11 11 (4-4,2-2) Ibrutinib gtable[layout]
## 12 12 (4-4,3-3) Ibrutinib gtable[layout]
## 
## [[2]]
## TableGrob (4 x 3) "Motolimod": 12 grobs
##     z     cells      name           grob
## 1   1 (1-1,1-1) Motolimod gtable[layout]
## 2   2 (1-1,2-2) Motolimod gtable[layout]
## 3   3 (1-1,3-3) Motolimod gtable[layout]
## 4   4 (2-2,1-1) Motolimod gtable[layout]
## 5   5 (2-2,2-2) Motolimod gtable[layout]
## 6   6 (2-2,3-3) Motolimod gtable[layout]
## 7   7 (3-3,1-1) Motolimod gtable[layout]
## 8   8 (3-3,2-2) Motolimod gtable[layout]
## 9   9 (3-3,3-3) Motolimod gtable[layout]
## 10 10 (4-4,1-1) Motolimod gtable[layout]
## 11 11 (4-4,2-2) Motolimod gtable[layout]
## 12 12 (4-4,3-3) Motolimod gtable[layout]
combiTab %>% filter(patientID == "P0033", drug2 == "QVD-Oph")
##  [1] wellID        rowID         colID         value         name         
##  [6] drug1         drug2         conc2         drug3         conc3        
## [11] concentration concIndex     obsViab       diagnosis     sampleID     
## [16] patientID     viab1         viab3         expViab       CI           
## [21] meanViab2    
## <0 rows> (or 0-length row.names)

Output per patient dose response curves - T-PLL (ibrutinib, motolimod)

drug.dfman <- drugList[["CombiManual"]]

drug.dfman <- drug.dfman %>% 
  separate(name, into = c("drug1", "drug2", "drug3"), sep = "[|]", remove = FALSE) %>%
  separate(drug2, into = c("drug2", "conc2"), sep = " ") %>%
  separate(drug3, into = c("drug3", "conc3"), sep = " ") %>%
  mutate(conc2 = str_replace(string = conc2, pattern = "ug/ml", replacement = ""), 
         conc2 = as.numeric(conc2))

drug.dfman$drug2 %>% unique()
## [1] NA           "CD95L"      "Adalimumab"
drug1Name <- c("Birinapant")
drug2List <- drug.df %>% filter(!drug1 %in% c("DMSO", "PBS", "empty"), !is.na(drug2), is.na(drug3)) %>% pull(drug2) %>%
  unique()

drug2List <- c("Adalimumab", "CD95L")

## Calculate the combination index for every drug
combiTab <- lapply(drug2List, function(x) {
  lapply(unique(drug.dfman$patientID), function(y) {
    drug.df <- drug.dfman %>% filter(patientID == y)
    
    message(x)
    # dose-response of drug1 (single agent)
    viabTab1 <- drug.df %>% filter(drug1 == drug1Name, is.na(drug2)) %>% dplyr::rename(viab1 = normVal)

    # dose-response of drug1 in combination with drug2
    viabTab2 <-  drug.df %>% filter(drug1 == drug1Name, drug2 %in% x, is.na(drug3)) %>% dplyr::rename(viab2 = normVal)

    # get the viability table of the drug with only 1 concentration (base drug)
    viabSingle <- drug.df %>% filter(drug1 == x, concentration == unique(viabTab2$conc2), is.na(drug2)) %>%
      distinct(patientID, normVal) %>% dplyr::rename(viab3 = normVal)

    testTab <- viabTab2 %>%
          left_join(select(viabTab1, patientID, viab1, concentration)) %>%
          left_join(viabSingle) %>%
          mutate(expViab = viab1*viab3) %>%
          mutate(CI = viab2/expViab, meanViab2 = mean(viab2)) %>%
          filter(!is.na(viab1), !is.na(viab2),!is.na(viab3)) %>%
    dplyr::rename(obsViab = viab2)
  }) %>% bind_rows()
}) %>% bind_rows()

lapply(c("Adalimumab", "CD95L"), function(x) {
  concInter <- drug.df %>% filter(name == "Birinapant") %>% pull(concentration) %>% unique()
  
  pList <- lapply(unique(combiTab$patientID), function(y) {
    combiTab <- combiTab %>% filter(patientID == y) %>% 
      filter(drug2 %in% x)
    
    if(x == "Adalimumab")  {lim.plot <- 1.35} else {lim.plot <- 1.25} 
  
    p <- combiTab %>%
      dplyr::rename(Expected = expViab, Observed = obsViab, Single = viab1) %>%
      pivot_longer(cols = c(Expected, Observed, Single), names_to = "viab",
                   values_to = "values") %>%
      mutate(viab = factor(viab, levels = c("Observed", "Expected", "Single")),
             drug2 = factor(drug2, levels = x), 
             patientID = paste0(patientID, " (", diagnosis, ")")) %>%
      arrange(desc(viab)) %>%
      #filter(viab %in% c("Observed", "Expected")) %>%
      ggplot(aes(x = concentration, y = values, group = viab)) +
      geom_line(linewidth = 0.5, aes(colour = viab)) +
      geom_point(size = 3.5, shape = 21, aes(fill = viab)) +
      geom_hline(yintercept = 1, linetype = "dashed", colour = "black", linewidth = 0.75) +
      geom_hline(yintercept = unique(combiTab$viab3), linetype = "dashed", colour = "darkorange", linewidth = 0.75) +
      scale_x_log10(labels = concInter, breaks = concInter) +
      scale_y_continuous(limits = c(0, lim.plot), breaks = c(0, 0.5, 1.0, 1.5)) +
      xlab("Birinapant (µM)") + ylab("Viability") + #ggtitle(x) +
      facet_wrap(. ~ patientID) +
      theme_classic() +
      theme(text = element_text(size = 22.5),
            plot.title = element_text(hjust = 0.5, size = 22.5),
            panel.background = element_rect(fill = "transparent",colour = NA),
            plot.background = element_rect(fill = "transparent",colour = NA),
            legend.background = element_rect(fill='transparent'),
            strip.background = element_blank(),
            #strip.text.x = element_text(size = 12.5),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.title = element_blank(),
            axis.line = element_line(linewidth = 1),
            #axis.text.y = element_text(size = 12.5),
            #axis.text.x = element_text(angle = 45, hjust=1, vjust = 1),
            #legend.position = c(0.25, 0.25
            legend.position = "none"
      ) +
      #scale_color_manual(values = c("Single" = "black", "Observed" = "#FFA000", "Expected" = "#4DB6AC"))
      scale_color_manual(values = c("Observed" = "#F44336", "Expected" = "#1976D2", "Single" = "black")) +
      scale_fill_manual(values = c("Observed" = "#F44336", "Expected" = "#1976D2", "Single" = "black"))
      # p

    #p
  })
  ggsave(plot = grid.arrange(grobs = pList, name = x), filename = paste0(opt$plot, x, "_combi_curve.png"),
         height = 35, width = 30, units = "cm")
  
  grid.arrange(grobs = pList, name = x)
})

## [[1]]
## TableGrob (4 x 3) "Adalimumab": 10 grobs
##     z     cells       name           grob
## 1   1 (1-1,1-1) Adalimumab gtable[layout]
## 2   2 (1-1,2-2) Adalimumab gtable[layout]
## 3   3 (1-1,3-3) Adalimumab gtable[layout]
## 4   4 (2-2,1-1) Adalimumab gtable[layout]
## 5   5 (2-2,2-2) Adalimumab gtable[layout]
## 6   6 (2-2,3-3) Adalimumab gtable[layout]
## 7   7 (3-3,1-1) Adalimumab gtable[layout]
## 8   8 (3-3,2-2) Adalimumab gtable[layout]
## 9   9 (3-3,3-3) Adalimumab gtable[layout]
## 10 10 (4-4,1-1) Adalimumab gtable[layout]
## 
## [[2]]
## TableGrob (4 x 3) "CD95L": 10 grobs
##     z     cells  name           grob
## 1   1 (1-1,1-1) CD95L gtable[layout]
## 2   2 (1-1,2-2) CD95L gtable[layout]
## 3   3 (1-1,3-3) CD95L gtable[layout]
## 4   4 (2-2,1-1) CD95L gtable[layout]
## 5   5 (2-2,2-2) CD95L gtable[layout]
## 6   6 (2-2,3-3) CD95L gtable[layout]
## 7   7 (3-3,1-1) CD95L gtable[layout]
## 8   8 (3-3,2-2) CD95L gtable[layout]
## 9   9 (3-3,3-3) CD95L gtable[layout]
## 10 10 (4-4,1-1) CD95L gtable[layout]

Output session info

## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
##  [1] stats4    parallel  grid      stats     graphics  grDevices datasets 
##  [8] utils     methods   base     
## 
## other attached packages:
##  [1] knitr_1.48                  scater_1.30.1              
##  [3] scran_1.30.2                scuttle_1.12.0             
##  [5] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
##  [7] Biobase_2.62.0              GenomicRanges_1.54.1       
##  [9] GenomeInfoDb_1.38.8         IRanges_2.36.0             
## [11] S4Vectors_0.40.2            BiocGenerics_0.48.1        
## [13] MatrixGenerics_1.14.0       matrixStats_1.3.0          
## [15] edgeR_4.0.16                limma_3.58.1               
## [17] gridExtra_2.3               circlize_0.4.16            
## [19] ComplexHeatmap_2.18.0       DrugScreenExplorer_0.1.0   
## [21] jyluMisc_0.1.5              lubridate_1.9.3            
## [23] forcats_1.0.0               stringr_1.5.1              
## [25] dplyr_1.1.4                 purrr_1.0.2                
## [27] readr_2.1.5                 tidyr_1.3.1                
## [29] tibble_3.2.1                ggplot2_3.5.1              
## [31] tidyverse_2.0.0             readxl_1.4.3               
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.3.2             later_1.3.2              
##   [3] bitops_1.0-7              cellranger_1.1.0         
##   [5] lifecycle_1.0.4           Rdpack_2.6               
##   [7] rstatix_0.7.2             doParallel_1.0.17        
##   [9] lattice_0.21-9            MASS_7.3-60              
##  [11] backports_1.5.0           magrittr_2.0.3           
##  [13] sass_0.4.9                rmarkdown_2.27           
##  [15] jquerylib_0.1.4           yaml_2.3.9               
##  [17] plotrix_3.8-4             metapod_1.10.1           
##  [19] httpuv_1.6.15             cowplot_1.1.3            
##  [21] RColorBrewer_1.1-3        multcomp_1.4-25          
##  [23] abind_1.4-5               zlibbioc_1.48.2          
##  [25] RCurl_1.98-1.14           TH.data_1.1-2            
##  [27] sandwich_3.1-0            GenomeInfoDbData_1.2.11  
##  [29] KMsurv_0.1-5              ggrepel_0.9.5            
##  [31] irlba_2.3.5.1             dqrng_0.4.1              
##  [33] commonmark_1.9.1          DelayedMatrixStats_1.24.0
##  [35] codetools_0.2-19          DelayedArray_0.28.0      
##  [37] xml2_1.3.6                DT_0.33                  
##  [39] tidyselect_1.2.1          shape_1.4.6.1            
##  [41] farver_2.1.2              viridis_0.6.5            
##  [43] ScaledMatrix_1.10.0       jsonlite_1.8.8           
##  [45] BiocNeighbors_1.20.2      GetoptLong_1.0.5         
##  [47] dr4pl_2.0.0               survival_3.5-7           
##  [49] iterators_1.0.14          systemfonts_1.1.0        
##  [51] foreach_1.5.2             tools_4.3.2              
##  [53] ragg_1.3.2                Rcpp_1.0.12              
##  [55] glue_1.7.0                SparseArray_1.2.4        
##  [57] mgcv_1.9-0                xfun_0.45                
##  [59] piano_2.18.0              shinydashboard_0.7.2     
##  [61] withr_3.0.0               fastmap_1.2.0            
##  [63] bluster_1.12.0            fansi_1.0.6              
##  [65] shinyjs_2.1.0             rsvd_1.0.5               
##  [67] relations_0.6-13          caTools_1.18.2           
##  [69] digest_0.6.36             timechange_0.3.0         
##  [71] R6_2.5.1                  mime_0.12                
##  [73] textshaping_0.4.0         colorspace_2.1-0         
##  [75] Cairo_1.6-2               gtools_3.9.5             
##  [77] tensor_1.5                markdown_1.13            
##  [79] utf8_1.2.4                generics_0.1.3           
##  [81] renv_1.0.7                data.table_1.15.4        
##  [83] htmlwidgets_1.6.4         S4Arrays_1.2.1           
##  [85] sets_1.0-25               pkgconfig_2.0.3          
##  [87] gtable_0.3.5              XVector_0.42.0           
##  [89] survMisc_0.5.6            htmltools_0.5.8.1        
##  [91] carData_3.0-5             fgsea_1.28.0             
##  [93] clue_0.3-65               scales_1.3.0             
##  [95] png_0.1-8                 km.ci_0.5-6              
##  [97] rstudioapi_0.16.0         tzdb_0.4.0               
##  [99] rjson_0.2.21              nlme_3.1-163             
## [101] visNetwork_2.1.2          cachem_1.1.0             
## [103] zoo_1.8-12                GlobalOptions_0.1.2      
## [105] KernSmooth_2.23-22        vipor_0.4.7              
## [107] pillar_1.9.0              vctrs_0.6.5              
## [109] maxstat_0.7-25            gplots_3.1.3.1           
## [111] slam_0.1-50               promises_1.3.0           
## [113] ggpubr_0.6.0              BiocSingular_1.18.0      
## [115] car_3.1-2                 beachmat_2.18.1          
## [117] xtable_1.8-4              cluster_2.1.4            
## [119] beeswarm_0.4.0            evaluate_0.24.0          
## [121] mvtnorm_1.2-5             cli_3.6.3                
## [123] locfit_1.5-9.10           compiler_4.3.2           
## [125] rlang_1.1.4               crayon_1.5.3             
## [127] ggsignif_0.6.4            labeling_0.4.3           
## [129] survminer_0.4.9           marray_1.80.0            
## [131] ggbeeswarm_0.7.2          stringi_1.8.4            
## [133] viridisLite_0.4.2         BiocParallel_1.36.0      
## [135] munsell_0.5.1             Matrix_1.6-1.1           
## [137] hms_1.1.3                 sparseMatrixStats_1.14.0 
## [139] statmod_1.5.0             shiny_1.8.1.1            
## [141] highr_0.11                rbibutils_2.2.16         
## [143] gridtext_0.1.5            drc_3.0-1                
## [145] exactRankTests_0.8-35     igraph_2.0.3             
## [147] broom_1.0.6               bslib_0.7.0              
## [149] fastmatch_1.1-4